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We make a detailed study of matter density perturbations in both metric and Palatini formalisms. 
Considering general theories whose Lagrangian density is a general function, f{R), of the Ricci scalar 
R, we derive the equation of matter density perturbations in each case, in a number of gauges, 
including comoving, longitudinal and uniform density gauges. We show that for viable f{R) models 
that satisfy cosmological and local gravity constraints (LGC), matter perturbation equations derived 
under a sub-horizon approximation are valid even for super-Hubble scales provided the oscillating 
mode (scalaron) does not dominate over the matter-induced mode. Such approximate equations are 
especially reliable in the Palatini formalism because of the absence of scalarons. 

Using these equations we make a comparative study of the behaviour of matter density pertur- 
bations as well as gravitational potentials for a number of classes of f{R) theories. In the metric 
formalism the quantity m = Rf,RR/f,R that characterises the deviation from the ACDM model is 
constrained to be very small during a matter era in order to ensure compatibility with LGC, but 
the models in which m grows to the order of 10~^ around the present epoch can be allowed. These 
models also suffer from an additional fine tuning due to the presence of scalaron oscillating modes 
which are absent in the Palatini case. 

In Palatini formalism LGC and background cosmological constraints provide only weak bounds on 
|m| by constraining it to be smaller than ~ 0.1. This is in contrast to matter density perturbations 
which, on galactic scales, place far more stringent constraints on the present deviation parameter 
m of the order of |m| < 10~^-10~'*. This is due to the peculiar evolution of matter perturbations 
in the Palatini case which exhibits a rapid growth or a damped oscillation depending on the sign of 
m. 



I. INTRODUCTION 

Recent high-precision observations by the Wilkinson Microwave Anisotropy Probe (WMAP) [l| together 
with high redshift supernovae surveys 0, observations of large scale structure Q and baryon acoustic 
oscillations JH have provided strong evidence that the Universe is at present undergoing a phase of accelerated 
expansion [a Such an accelerating phase poses a serious problem for cosmology since it is difficult to 
explain within the usual general relativistic framework. 

Phenomenologically, the simplest way to generate such an accelerating phase is through the addition of a 
cosmological constant to the Einstein's field equations. At a more fundamental level, however, the problem 
is how to account for such a constant within a candidate theory of quantum gravity. This has motivated a 
large number of alternative models (see Ref . @ for a recent review) . These models can mostly be divided 
into two broad groups: those involving an exotic matter source and those modifying the gravitational sector 
of the theory. An important subset of the latter are the so called generalised gravity theories, based on 
non-linear lagrangians of the form /(i?), where / is a general differentiable function of the Ricci scalar R. 
Such modifications to the (linear) Einstein-Hilbert action typically arise in effective actions derived from 
string/M-theory 

There are two approaches that may be taken in studying generalised f{R) theories of gravity, depending 
upon the choice of the independent fields with respect to which the action is varied. In the first ('met- 
ric' approach) only variations with respect to the metric are considered, whereas in the second ('Palatini' 
approach) the action is varied with respect to both the metric and connection. Both approaches result 
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in identical field equations for the case of linear Einstein-Hilbert action. In the more general cases, with 
nonlinear Lagrangians f{R), however, the metric approach results in fourth-order field equations, whereas 
the Palatini variation generates a second-order system. In this paper we consider both approaches. 

Recently a number of /(i?) models have been proposed in order to explain the late-time acceleration of 
the Universe. Given the large number of such models that have been (or can potentially be) considered, 
an urgent task at present is to devise stringent tests in order to reduce the viable range of candidates. In 
addition to constraints obtained by demanding the stability of the models, we require the models to be 
consistent with constraints from background cosmological dynamics [ill , [l^ , as well as from local gravity 
constraints (LGC) [H, [H, [H [3 . 

Generally, models put forward to explain the late-time acceleration involve infra-red corrections to Hilbert 
action with negative powers of the Ricci scalar R. It is interesting to note that the viability of such models 
can depend upon which formalism is used in order to derive them. As an example consider the model 
f{R) = R - with n > [U [3 (see also Ref^_pi). In the metric approach this model is 

problematic because of the absence of the matter era [ll|, |20|| as well as the instability of perturbations 
associated with negative f b.r [HI [131, where comma denotes differentiation with respect to R. Moreover it 
does not satisfy the LGC [l3| . 

On the other hand, in the Palatini approach the background cosmological dynamics successfully generates 
a standard matter era followed by a late-time acceleration [1^. One can also realise a correct Newtonian 
limit in the regime where R is much larger than [24] . For the study of cosmological and local gravitational 
constraints for a range of f(R) models, see Refs. in the metric case and Refs. [1^ in the Palatini case. 
See also Refs. [13, [H, [H, Ho] for f{R) theories that satisfy both cosmological and local gravity constraints 
in the metric approach. 

Despite the importance of above constraints in limiting the range of viable models, it is well known 
that the study of the homogeneous dynamics is not sufficient on its own to determine the nature of the 
underlying theory uniquely, given that f{R) theories can always be expressed in terms of a conventional 
relativistic cosmology sourced by a perfect fluid with an effective equation of state. An important way 
to break this degeneracy is by considering the evolution of density perturbations about the Fricdmann- 
Lemaitre-Robertson- Walker background. The study of perturbations provides an important tool in order to 
break this degeneracy as well as allowing more stringent constraints to be placed on the parameters of the 
models. In fact a number of authors recently studied the evolution of density perturbations for some /(i?) 
models to put constraints on model parameters in metric [H, [3l|, [H, [111 and Palatini [3l,[3l,[3^ formalisms. 

Here, with this aim in mind, we study the evolution of density perturbations and resulting observational 
consequences in both metric and Palatini formalisms. We first write the equations without specifying any 
gauges and then write them in a number of gauges (including comoving, longitudinal and uniform density 
gauges) from which we derive the corresponding approximate perturbation equations under sub-horizon type 
approximations. In the metric approach where the oscillating mode (referred to as the scalaron, see below) 
[37| is present, this approximation can be invalid if the scalaron is overproduced in the early Universe. 
However, as long as the scalaron is sub-dominant relative to a matter induced mode, we shall show that 
approximate perturbation equations can be valid even for super-Hubble modes in the models that satisfy 
LGC. The approximation is especially reliable in the Palatini case because of the absence of scalarons. 
The simplicity of the equations derived facilitate the estimation of the growth rate of perturbations both 
analytically and numerically. 

Using these equations we make a comparative study of the behaviour of matter density perturbations in 
both formalisms, for a number of classes of theories satisf ying LGC as well as background constraints. These 
include viable theories recently proposed in the literature [27l [28l . [29j . An important quantity to characterise 
the growth rate of matter perturbations is s = S'^/Sm, where Sm is the density contrast on orthogonal 
comoving hypersurfaces and a prime represents a derivative with respect to the number of e-foldings. In the 
standard general relativistic (GR) case s — 1. 

In the metric formalism, the growth rate in the non-standard regime where ^ = kF' / [a^R) to ^ 1 (a is the 
scale factor and k is the comoving wave number) is given by s = 1.186. While this is within the current 
observational bound s ^ 1.5, the difference between spectral indices of the matter power spectrum and the 
CMB spectrum can provide stronger constraints. In the Palatini case, when to > 0, the growth rate s grows 
exponentially once the Universe enters the regime ^ ^ 1. Thus the quantity to can be severely constrained 
using the information from the bounds on s. In this paper we shall obtain constraints on the present value 
of the deviation parameter to as well as model parameters for a number of /(i?) theories that include most 
essential features required by viable models. 

The plan of the paper is as follows. In Section |TT] we give a brief review of constraints for f{R) theories 
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in both metric and Palatini formalisms provided by background cosmological evolution and local gravity 
constraints. In Section [IIII we give a brief account of metric perturbations and matter density perturbations 
which are applicable to both metric and Palatini formalisms. In Sections IIVI and fVl we derive the evolu- 
tion equations for matter density perturbations in metric and Palatini formalisms respectively and in each 
case discuss the constraints they provide by considering a number of theories that are viable according to 
background cosmological and local gravity constraints. Finally we conclude in Section IVll 

In what follows we shall use units such that SttG — 1, but we restore bare gravitational constant G when 
it is needed. 



II. LOCAL GRAVITY CONSTRAINTS AND VIABLE BACKGROUND COSMOLOGICAL 

EVOLUTION 



The classes of f{R) general gravity theories we consider have actions of the form: 

„ 1 



d'^Xy^ f{R) + Srn{gfi,y,1pm) 



(1) 



where / is a general differentiable function of the Ricci scalar R and Sm corresponds to a matter action, 
which depends upon the metric and matter fields ipm- The Ricci scalar R is defined by i? = g^'^Rfj.u, 
where the Ricci tensor _R^j/ is 



(2) 



In the case of the metric formalism, the connections P^^ are the usual metric connections defined in terms 
of the metric tensor 9^,^. In this case the field equations are obtained by varying the action ([T]) with respect 
to g^i, to give 

1 



F{R)R^,{g) - -fiR)g^, - V^V,F(i?) + g^,nF{R) = T^, , 



(3) 

SUFiR) + F{R)R-2f{R) , (4) 

where F — df /dR. Note that Eq. ^ corresponds to a trace of Eq. ([3]). The third and fourth terms on the 
left hand side of Eq. ^ come from the variation of -R^^y with respect to g^^u. 

In the case of the Palatini formalism P^^ and g^^ are treated as independent variables. Varying the action 
dJ) with respect to g^j, gives 



F{R)R^,{r) ~ = T,, , 



(5) 



where (P) is the Ricci tensor corresponding to the connections P^^ , which is in general different from the 
Ricci tensor corresponding to the metric connections i?^^ (g) . Taking the trace of this equation, we obtain 

F{K)R-2f[R)^T, (6) 

where R{T) = g^'^R^^iT) is directly related to T and is different from the Ricci scalar R{g) — g^'^ R^„{g) in 
the metric case. Taking the variation with respect to the connection, and using Eq. ([5]), we find 



~ -i^gt^vRig) ^ ^ ^ Sm-^ 



d^Fd,F--g^,iVFf 



•(7) 



It is useful to express the generalised theories based on action ([T]) in terms of generalised scalar-tensor 
theories of Brans-Dicke type with a scalar field cj) and the corresponding potential ^(0) thus [38l |: 



-g 



+ Smigfiiy, V'm) : 



(8) 



where ojbd is the Brans-Dicke parameter. Note that the original Brans-Dicke theory [3^ corresponds to 
V{(f) = 0. To see the correspondence of the action (|8]) with ([1]), we recall that the variation of ([8]) with 
respect to gf^i, and (j) leads to the following equations 



Ri^uig) - i^g^cuRig) 

(3-f 2u;BD)n0 + 2y(0) 



-T . 



1 

d^,(l)d^(j>- -g^A^4>y 



,(9) 
(10) 
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Now consider the following correspondence 

^^F{R), V{^) = Ri(t>)F^f{Ri4>)), (11) 

where R = R{g) in the metric case and R — R{T) in the Palatini case. Comparing Eqs. ©-(HUl) with 
Eqs. ©-(HI), it can readily be seen that the f{R) theories in the metric formalism correspond to the above 
generalised Brans-Dicke-type theories with wbd = 0. Similarly the f{R) theories in the Palatini formalism 
correspond to the Brans-Dicke-type theories with wbd = —3/2. 

The Newtonian effective gravitational constant in the Brans-Dicke theory ^ can be derived under a 
weak- field approximation by considering a spherically symmetric body with a mass Mq, constant density p 
and a radius r© and a vanishing density (p = 0) outside the body. Decomposing the field (j) into background 
and perturbation parts (0 = 00 + ^4') ^nd using a linear perturbation theory in the Minkowski background 
with a perturbation /i^j^, the effective gravitational constant is given by [fSj]: 

^Newton = h + , for Mtq « f . (f2) 



3 -I- 2wbd 

Here ^ is a distance from the center of the body and the scalar field mass squared is defined by 

where is a local field in Minkowski spacetime. We should emphasise here that the expression (If 2p is only 
valid subject to Mr© <C f ff the scalar field mass M is large, it can happen that the condition for the 

linear perturbation theory {dip <C 4>o) becomes invalid. Moreover, this validity depends on the distribution 
of scalar field mass inside and outside the body. When the mass in the region £ < r© is much larger than the 
corresponding mass in the region £ > r©, a thin-shell can be formed inside the body to satisfy local gravity 
constraints through a chameleon mechanism 40]. The formation of the thin-shell occurs in a non- linear 
region in which the above linear result p2)) ceases to be valid f4]. 

An important point to note here is that the Palatini case, corresponding to ciJbd = —3/2, is rather special 
in a number of fundamental ways. For example, the cj) field kinetic term in Eq. ()f 0|) vanishes in this case, 
whereas it is non-zero in the metric case with wbd = 0. As we shall see below this has the important 
consequence that the oscillatory scalaron mode is absent in /(i?) theories based on the Palatini formalism, 
whereas it is present in all other models with cijbd —3/2, including f{R) theories based on the metric 
formalism. 

Also the scalar field mass M blows up for finite potential-dependent terms in the parenthesis of Eq. 
as wbd approaches —3/2. For the theories with wbd —3/2, the scalar field mediates a fifth force with an 
interaction range f/M because the field has a kinetic term in Eq. (fTU|) . In the Palatini formalism the field 
lacks a kinetic term in Eq. pop , which implies that the notion of the usual interaction range determined by 
mass M does not hold. Thus the Palatini case should be treated separately compared to other theories. 

In the usual Brans-Dicke theory with V{(t>) — and ljbd —3/2, the mass M vanishes because the 
field (f) propagates freely. Then the Yukawa-correction term e~^^^ in Eq. (|12l) becomes one, in which case 
the Brans-Dicke parameter wbd is constrained by local gravity experiments to be larger than 40000 (4l| . 
This constraint does not, however, necessarily apply to f{R) gravity theories, because the presence of the 
potential V{(j)) can make such theories compatible with local gravity constraints under certain conditions 
(which we shall see later). 

Before proceeding to consider matter density perturbations and the constraints arising from them, we shall 
first briefly review the cosmological and local gravity constraints for /(i?) theories in both formalisms. The 
constraints from background dynamics have been studied in both metric [if], [H, HO] and Palatini H, 34 1 

m 



formalisms. To discuss these constraints, it is useful to introduce the following dimensionless parameters 



Rf.RR Rf.B 

m = —J , r = — . (14) 

J,R J 

In the metric case insight into the cosmological dynamics can be obtained by considering trajectories in 
the (r, to) plane. In the Palatini case it is more convenient to use other dimensionless variables for the 
background dynamics [l^, but the deviation parameter to is still important in order to characterise the 
deviation from the ACDM model. 
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A. Metric formalism 



We first review the cosmological viability in the metric case and then proceed to consider the LGC. In the 
(r, to) plane the matter point corresponds to Pm '■ ir^m) « (—1,0). In order to have a saddle matter era 
followed by a late-time acceleration, we require the following conditions [ll| 

m>0, -1<^<0, at (r,m) « (-1,0) . (15) 
dr 

De-Sitter fixed points Pa lie on the line r — —2. They are stable provided that 

r = -2, 0<m<l. (16) 

As long as conditions ([T5|) and (fTB|) are satisfied, we can realize a successful cosmological sequence (radiation, 
matter, and de-Sitter epochs). 

There are a number of models in the recent literature that satisfy the above cosmological constraints. 
Examples are (i) f{R) = a{R'' - with c > 1, 6c « 1 132*1, and (ii) f{R) = R - aR<^ with a > and 
< /9 < 1 [l2|. For these models, the parameters to and r satisfy the relation to = C{—r — 1), where C 
is a positive constant in the neighbourhood of r = — 1. Using observational constraints on the background 
dynamics from SN la and the sound horizon of CMB, the parameter to has been shown to be constrained 
to be TO < 0(0.1) 

If the information from LGC is also included the constraints on the model parameters become very strong. 
The usual procedure to determine the local gravity constraints for /(i?) theories is to consider their Brans- 
Dicke representations (jlip and expand the equations of motion around a background Minkowski metric [T3| . 
Since the metric formalism corresponds to ujbb = Oi the scalar field mass defined in Eq. (jl3p is given by 

.r2 1 f f,R 



where we have used the relation PT|) . If < the Yukawa correction e~^^^ is replaced by an oscillating 
function cos(|M|^), but this case is excluded by the experimental requirement that 7 « 1. Hence the mass 
squared Af^ is required to be positive. 

Clearly we require the mass M to be heavy for consistency with local gravity experiments. In that case, 
however, the effective gravitational constant (|12p obtained under the linear approximation ceases to be valid. 
As was already mentioned above, a thin-shell begins to form through a chameleon effect in this non-linear 
regime. To consider this chameleon effect in f{R) gravity, it is convenient to introduce a new conformally 
related metric and a scalar field [4^ : 

Qtj.i' = (l>9iiu , if = ^/^\Ti(t>. (18) 
Then the action llj in the Einstein frame takes the form: 

'1^ I 



S = / d^XA/^ 



^R--{V^f-U{^) 



+ ^m(gM-e^'^'^,V„0, (19) 



where the coupling (3 in /(i?) models and the potential U are given by 



In a spherically symmetric setting with an energy density p (= — T), the field ip satisfies the following 
equation Q 

dV^2d^^dKff .21) 

where (. is the distance from the center of symmetry in the Einstein frame and 

Feff(^) = [/(^)-He^V- (22) 
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Here we have introduced an energy density which is conserved in the Einstein frame |40| ]. 

As an example, let us consider a spherically symmetric body that has an energy density p* = p\ inside 
the body < f© = e^^'^r©). Let the energy density outside the body (l> f©) be given by p* = P*b, which 
is much smaller than p\. Then the effective potential (|22p has two minima at = (pA and p = ps satisfying 
the relations 



U^iipA) + Pe^'^^p*A^O, 



(23) 
(24) 



The effective mass at the potential minima are given by m\ = V^g{pA) S> = V"g{pB), which comes 
from the condition P*a^ P*b- 

Imposing appropriate boundary conditions at £ = and £ = f©, the solution to Eq. (|2ip can be approxi- 
mated by [i3, liO, Si] 



p{£) 

where M© = 47rr|pA/3 = 47rf|p^/3, 



47r 



/3cff = Ap— 



Af© 

^0 



6/3$ 







(25) 



(26) 



and $© = GMq/tq. A thin-shell is developed under the condition Af©/f© <C 1. In this case the effective 
coupling I /3eff I becomes much smaller than \[3\ so that the models can be consistent with local gravity 
constraints |40 |. 

Models that can satisfy these conditions have recently been proposed by (i) Hu & Sawicki [23| and (ii) 
Starobinsky jH]: 



(i) f{R) = R-\R, 



{R/Rcf 



{R/Rc 



\2n 



1 ' 



(ii) f{R) = R-XR, 



1 - 



(27) 
(28) 



where n, X and Rc are positive constants. In both models the cosmological constant disappears in a flat 
spacetime, i.e., /(i? = 0) = 0. Note that Rc is roughly of the order of the present cosmological Ricci scalar 
i?o for n = 0(1) and A = 0(1). In high curvature regimes R^ Rc these models behave as 



/(i?) ^R-XRc 



R 



(29) 



with 



m ' 



C{-r-l) 



2n + l 



(30) 



where C is a positive constant, and m and r are defined in Eq. Thus they are very close to the ACDM 

model with suppressed values of m during matter and radiation eras (r ~ — !)• 

In the regimes R ^ Rc one can show that the term \pB — Pa\ in Eq. (|26p is of the order of m{RB) for 
n = 0(1), where Rb is the Ricci scalar in the neighbourhood of ipB (which is generally much larger than Rc 
in an environment where local tests of gravity are carried out). Hence the thin-shell is developed under the 
condition 



(31) 



This can be regarded as a criterion for the compatibility with local gravity constraints. In the case of the 
earth, the condition ([3T|) corresponds to m{RB) ^ $© ^ 10^^. Since <&© ^ 1 in most of local gravity 
experiments, the parameter m is constrained to be much smaller than the order of unity in the region where 
the Ricci scalar Rb is much larger than the present cosmological one {Rq ^ Rc). 
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Cosmologically the condition (PT|) implies that viable models need to be very close to the ACDM model in 
radiation and matter dominated epochs {R ^ Rq). However, deviations from the ACDM model are allowed 
around the present acceleration epoch (i? ^ i?o)- Thus in viable models the parameter m is negligibly small 
during the radiation and matter eras, but continues to grow by the present epoch. 

For the theories of the type ((29)) the corresponding Brans-Dicke field (j>, the potential V{(f>) and the mass 
squared are given by 



1 - 2n\{RjR) 



2n+l 



1/(0) ~ XR, 



1 - (2n+ 1) 



2n\ 



Rr 



3(2n+ 1 



■(2nA)^^(l 



(32) 
(33) 

(34) 



which in the limit R/Rc oo become 0^1, V{(j)) XRc and oo respectively. In these regimes 

the field is stuck around (/) — 1 because of the presence of a p-dependent term. When R decreases to the 
order of Rc, the field begins to evolve along the potential V{(j>) with a lighter mass M which is not very 
much different from Rc- Thus in the Brans-Dicke description, the departure from the point 0=1 amounts 
to deviation from the ACDM model. 

The models ([27)1 and (p8)) are constructed to satisfy the stability conditions 



/,i?>0, f,RB.>0, for R>Ri{>0) 



(35) 



where i?i is a Ricci scalar at the late-time de-Sitter point. The first condition is required to avoid repulsive 
gravity, whereas the second ensures the absence of tachyons or ghosts. The second condition is also required 
for the consistency with LGC (as was shown above) as well as to ensure the stability of density perturbations 
[23 , [3H (as we shall see below) . We also note that the requirements (I5T|) and ([55)1 are perfectly consistent 
with the condition < m{R) <C 1 derived in Ref. [ll| for the existence of the matter era. 

To summarise, the conditions ([3T|) and (|35|) . together with the existence of the de-Sitter point ([TH]) . are 
required for the viability of f{R) models in the metric formalism. The condition for the existence of the 
saddle matter era given in Eq. is automatically satisfied under the requirements (Plj) and 



B. Palatini formalism 



Let us next consider cosmological and local gravity constraints for /(i?) theories in the Palatini formalism. 

We first discuss the conditions for the cosmological viability of f{R) models at the background level. In 
Ref. [235 it was shown that radiation {Pr), matter (Pm) and de-Sitter (Pa) points exist irrespective of the 
forms of f{R) provided that the function 

''^''^ - -\FR^f)iFnR~F) ' 

is well-behaved (i.e., it does not show discontinuous or divergent behaviour). Note that effective equations 
of state corresponding to points Pr, Pm and Pa are given by Wcs — 1/3,0,-1, respectively. The de-Sitter 
point Pa corresponds to FR — 2f — 0, i.e., 

r = -2, (37) 

and C{R) = 0. This de-Sitter solution exists on the same line as in the metric case. Since the eigenvalues 
of the Jacobian matrix for perturbations around the point Pd are (Ai, A2) = (—3 — C{R), —4 — C{R)), the 
de-Sitter point on the line r = — 2 is always a stable attractor. This situation is different from the metric 
case in which the stability of the de-Sitter point requires the additional condition < m{r = —2) < 1. 

The stability of radiation and matter points, on the other hand, depends upon the particular f{R) models 
chosen. The ei genv alues for perturbations are given by P^: (Ai,A2) = {i + C{R),l) and Pm- (Ai,A2) = 
(3 -I- C{R), —1) |23|. Thus the models with C{R) > —3 give rise to the unstable node for Pr and the saddle 
point for Pm- Thus models satisfying the condition C{R) > —3 lead to a sequence of radiation, matter and 
de-Sitter epochs. As an example, let us consider the following model [l3, [3 

2(n+l) 

fiR) = R-^^^^ (38) 
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where /i and n are constants. In this case one has C{R) = 3n in the regime R^~^^ ^ which means 

that a successful background trajectory is reaUsed for n > —1. Note that a stable de-Sitter solution exists 
with RI'^'^ = (2 + n)/i^("+^) and C{R) = 0. Obviously the constraints for a successful trajectory, at least at 
the background level, are not so severe compared to the metric formalism. One does not even require the 
condition to > for the existence of a viable matter era. 

Let us next discuss LGC in the Palatini case. Considering the Brans-Dicke picture of f{R) theories, 
the theories in the Palatini formalism correspond to ujbd = —3/2. Now since the usual description of 
the interaction range determined by the inverse of mass M can not be applied in this case, one needs to 
proceed in a different wav|T3|. From Eq. ^ the field (j) = depends upon the value of the trace T, 

i.e., (p = We expand the field around the vacuum: (j){T) = 00 + {dr^T + ■ ■ ■, where (po = (j){T = 0) 

and T « — p[l — 0{v'^ /c^)]. Note that we use the non-relativistic approximation under which the velocity v 
of matter is much smaller than the speed of light c. Carrying out a post-Newtonian expansion around the 
Minkowski vacuum {gfn, — rj^^^ + hf^^) in the solar system, the solutions for the second-order perturbation 
equations are given by 



(39) 
(40) 



.(2) _ 2GN-t°"M0 Vo 


"oo — 


i 6(/)o ■ 


"-v - 


■27GN™to„^j^ T/o 


i 600 



where Vq 



V{(j)o)- Here the effective gravitational constant and the post- Newtonian parameter are 

Mq - Mv 



G Newton 
off 



00 



My 

Me 



1 



Mq+Mv 



where Mq and My are given by 



d^xp(i,x)- 



Mv 



d^x 



Vo 

00 



V 



(41) 



(42) 



and p is the energy density of the sun. 

To ensure LGC, three conditions need to be satisfied [if 



. (i) \Mv\ « iMol, 
. (ii) \Vae/M « 1, 



• (iii) The contribution of the term log (0/0o) is negligible. 

The first condition comes from the experimental requirement 7 « 1. Since it is not easy to interpret this 
requirement directly, we shall elucidate this by considering a specific f{R) model later. Concerning condition 
(ii), setting T = in Eq. ^ and using pT|) to obtain Vq = f{Ro), this translates into 



fiRo) 



r < 1. 



(43) 



When the deviation from the ACDM model is not large, the term /(i?o)//,i?o(^o) is of the order of the 
present cosmological Ricci scalar Rq ^ Hq. Hence on the scales of the solar system this condition is well 
satisfied. 

Regarding condition (iii), the presence of the term log (0/0o) in Eqs. and (|40p leads to an additional 
acceleration of particles that should be small in order to be consistent with experiments. From the validity 
of classical Euler equations, the condition (iii) translates to [Tsj 



Psd^/dT 



« 1, 



(44) 



where ps is an energy density of the local structure, 
strong dependence on T. Using the relations T = 2V - 
translates to 



This means that the field 0(T) should not have a 
and 0T^00 — V4, = f,B./ f,RR — R, this condition 



PsII.. 



Rs 



Irs/1 



RsRs 



R 



< 1. 



(45) 
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It is sometimes useful to rewrite this condition in terms of the variable m thus: 

- 1 



m{Rs 



When \m{Rs)\ ^ 1, this is well satisfied since both /j^^ and Ps/Rs are of the order of unity. Note that this 
constraint is not so restrictive compared to the metric formalism. This can be understood by recalling that 
in the Palatini case the field is non-dynamical without an interaction range. In the metric formalism one 
needs a large scalaron mass M to satisfy the thin-shell condition, which leads to a very small value of m{Rs) 
satisfying Eq. (pij) . We also note that in the Palatini case the condition J.rr > is not required in order to 
satisfy LGC. 

As a concrete example, let us apply the above constraints to the theories given by Eq. ([38]) with n> —1. 
In order to give rise to a late-time acceleration, p, needs to be of the order of the present Hubble radius Hq. 
The field </) and the potential V{(f)) defined in Eq. (fTTj) are in this case given by 



0=l + n(^^j ' 

= + ^l'^{n + l)^,^[t-L^ . (48) 

Now in the de-Sitter case [vacuum (T = 0)] the solution i?o satisfies 

F{Rn)Rn~2fiRo) ^0, (49) 

which for the model ([38]) gives 

i?o = (n-H2)i/("+i)/x2, (50) 

, 2{n + l) n+1 2 

'^" = ^rT^' ^°^(n + 2)«/("+i)'^ • 

In settings where local gravity experiments are carried out, the parameter e = ^ Po/Ps is much 

smaller than unity. For example if we take the mean density ps — lO^^^g/cm'^ and use the typical values 
p^ ^ Hi ^ po = 10~^^ g/cm^ and Rg ^ ps, then e is of the order of 10~^^. 

When 71 > 0, then in the limit e 0, we have ^ ^ 1 and V{(t>) 0. Thus in the expression of Mq given 
in Eq. the term Vo/0o dominates over the term V/cf), thus giving 



and 



Mv « j d^xVo^ J d^x p^ , MqK j d^x ps 



(52) 



Now since p'^ ^ po <^ ps, then the condition (i) is well satisfied. 

When — 1 < n < 0, then as e — + one has </)—!■ 1 and the potential V becomes of the order V ~ 
p^ip^/R)" :s>Vo^ p^. This gives 

\Mv\ « J d^xp^p^R)" - J d^xpoipo/psT , (53) 

where Mq is the same as that in Eq. ([5^ . The ratio of the integrands in the expressions for My and Mq 
can be estimated to be (po/ps)"^^ ^ 1, which means that the condition \Mv\ ^ Mq is again satisfied. 
The parameter m in this case is given by 

(n + l)ne"+i , , 

Now since e is much smaller than 1 we obtain \m{Rs)\ ^ 1. Hence theories of type ([55]) with n > —1 can 
satisfy local gravity constraints. 

The above discussion shows that it is easier to satisfy the local gravity constraints in the Palatini case 
than in the metric case. In the latter case we also require the condition J^rr > to ensure that the scalaron 
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mass squared is positive. Moreover the requirement of the heavy mass M gives very small values for 
m{Rs), which imposes the condition that viable f{R) models need to be very close to the ACDM model 
during matter and radiation epochs. We also note that even though the condition \m{Rs)\ ^ 1 is also 
required in the Palatini case, the absolute values of |TO(i?s)| do not need to be vanishingly small. In fact 
even models (j38p with n > can satisfy the correct Newtonian limit, while they are excluded in the metric 
formalism because f.RR is negative in those cases. Thus in the Palatini formalism the models of the type 
f{R) = R — g{R) can be consistent with local gravity tests provided that the contribution of the term g{R) 
is not significant relative to the linear term. 

In subsequent sections we discuss the evolution of density perturbations for /(i?) theories in both metric 
and Palatini formalisms. We shall carry out a detailed analysis for a number of f{R) models that can satisfy 
both the cosmological and local gravity constraints and use the evolution of density perturbations to place 
constraints on the model parameters as well as their deviation from the ACDM model. 



III. MATTER PERTURBATIONS AND GAUGE ISSUES 



In this section we present the equations for matter perturbations that are applicable to f{R) theories 
in both metric and Palatini formalisms. As our background spacetime we shall consider a flat Friedmann- 
Lemaitre- Robertson- Walker (FLRW). The perturbed FLRW metric which includes linear scalar metric per- 
turbations a, b, If and E can be written in the form [43 | 

ds^ = -(1 + 2a)dt^ - 2a5,,dMa;* + a{t)^ [(1 + 2^)% + 2Ei,j] dx'dx^ , (55) 

where a{t) is a scale factor. We shall consider a pressure-less matter source. Now since the flow is irrotational 
for scalar perturbations, we can introduce a velocity potential V in terms of which the components of the 
energy momentum tensor of the pressure-less matter can be decomposed as 

Tq = -{Pni + Spm) , T" ^ Pm(y -b)^i = ~PmVm,i, (56) 

where Vm is related to the velocity potential V through [4^ : 

Vm^-iV^b). (57) 

Note that the definition of b in Ref. [i^ has an opposite sign to that used here. We recall that in both metric 
and Palatini formalisms the matter energy density pm satisfies the standard continuity equation 

Pm + iHp^ = , (58) 

where H = d/a is the Hubble parameter and a dot represents a derivative with respect to t. The matter 
perturbation can then be shown to satisfy the following equations of motion in the Fourier space [4^ |47| : 



Spm + SHSpm = P7n(^K~ 3Ha - -^""m^ ; (59) 
+ Hvra = -a , (60) 



a 



where k is a comoving wavenumber and 



fc2 

3iHa-^) + —X, x=a{b + aE). (61) 



Defining the following variables 



a«™ = -a(l/ - 6) , ,5=^, (62) 

Pm 



where w is a covariant velocity perturbation [48l|. Eqs. ((59)) and ([60)) can be written as 



a = V , (63) 
5^K~ 3Ha - —V . (64) 
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Now choosing a comoving hypersurface, the density perturbation can be expressed in a gauge-invariant 
way as [isj : 

5p,n^ 5p,n + af}yyi{V -b) . (65) 
We shall define the density contrast on comoving orthogonal hypersurfaces as 

<5„, = ^ + 3iJz;. (66) 

Pm 

Now since the right hand side of Eq. is gauge-invariant, 5m can be evaluated in any gauge. The evolution 
equation for 5m is given by 

5m + 2H5m + ^{a~x)^^B + QHB, (67) 

where B = Hv — (p. 

In the following we shall consider three different gauges to fix the gauge degree of freedom: 

• Comoving gauge: in which the 3- velocity and the scalar shift function vanish (i.e., w = 0). This 
implies that along with the 3-velocity the momentum vanishes as well. Thus the gauge-invariant 6m 
in this gauge becomes 



^ Spm 
Pm 



(68) 



• Longitudinal gauge: in which the shift vector b and the anisotropic potential E both vanish, resulting 
in X = 0. The gauge-invariant 5m in this gauge becomes 



(69) 



• Uniform density gauge: in which we have constant density hypersurfaces, i.e. 5pm = 0. The 
gauge-invariant 5m in this gauge becomes 



<5(f ^ = 3Hi 



In what follows we discuss the evolution of matter perturbations in the above three gauges, while noting 
that the physics does not depend upon the choice of gauges. Also while the above discussions hold in both 
metric and Palatini formalisms, the equation of matter perturbations is different in each case as we shall see 
in the following sections. 

IV. DENSITY PERTURBATIONS IN THE METRIC FORMALISM 

In the metric formalism the background equations are given by 

3Fij2 ^^(FR-f)- 3HF + pm , (71) 

-2FH = F - HF + Pm, (72) 

where R — 6{2H^ + H). Note that we only take into account a non-relativistic matter. In Fourier space the 
scalar metric perturbations satisfy the following equations of motion, in the so-called gauge-ready form [46^1 



'^ip + 3H(Ha -p) + -^Hx = ttt^ SH5F - { 3H + 3H^ - ^ ] 5F - 3HFa -Fk- 6pm , (73) 

zr \ a-' 
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Ha — if 



1 



2F 

X + Hx- a- LP 



5F - H5F -Fa + p^v 
1 



F 



{SF - Fx) , 



k + 2Hk 
1 

^ 2F 



-6H^ + ^]6F + SHSF + S5F - Fk - 3(2F + HF)a - 3Fa + 5p„ 



SF + 3HSF 



R 



SF = \sp„, + F{k + a) + {2F + 3HF)a - \f6R . 



(74) 
(75) 

(76) 
(77) 



In the following we derive the perturbation equations in the above gauges, in both exact forms as well as 
using a sub-horizon approximation. 



A. Comoving gauge (« = 0) 



We first derive the equation of matter perturbations in the comoving gauge {v = 0). When d = we have 
a = and Sm^ = k from Eqs. ([55]) and Hence from Eq. ([7S|) we find 



2H 



2F 



1 

2F 



SF + 3HSF + 3SF + Sp„ 



(78) 



whereas from Eq. ([77)) . the perturbation SF satisfies 



SF + 3HSF + - 4ff2 _ 2ij) SF = lsp„, + FS^^^ . 



3/, 



RR 



(79) 



The evolution of the matter perturbations Sm"^ can then be obtained by solving Eqs. ([78| and ((79l) numer- 
ically. For models that satisfy local gravity constraints the mass squared term defined in Eq. ()17p can be 

r2 ^ f,R 



well approximated by M 



a term which appears on the left hand side of Eq. ([79|) . We are mainly 



interested in the evolution of perturbations on sub- horizon scales, i.e. 



fc2 



:>{hMh\}. 



(80) 



We also recall that for the models that satisfy LGC the mass of the scalar field squared is much larger 
than R ~ \H\. Hence either k'^ /a? or A'P is dominant in the parenthesis on the left hand side of 

Eq. ((79|) . Let us first consider the case in which the time-derivative terms in SF are neglected, i.e.. 



-\SF\.AI^\SF\ ■> {\HSFl\SF\} . 



(81) 



The condition ([5T|) amounts to neglecting the term SF that leads to the oscillation of SF . This is the 
approximation used in scalar tensor models in Refs. [1, [4^ [13] ■ Later we explore the validity of such an 
approximation paying particular attention to the conditions that should be satisfied. 
Under the conditions dSD]) and dH]), Eq. dUD gives 



SR^ 



1 Sp,n + 'iFS : 

F l + iS, 



[v) 



where 



f.RR 

a' f.R 



a'^R 



-m . 



(82) 



(83) 
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Using the approximation ([5T|) in Eq. (|78p . we obtain 



4") +I2H+ ^ j Si:^ - 47rG^|^-°p„,5(„") ^ , (84) 



where the "cosmological" effective gravitational constant is given by 



= S f^^l ■ (85) 



F + 

Note that we have restored the bare gravitational constant G. 

Introducing a physical wavelength £ = a/k, the parameter ^ defined in Eq. (|83p can be written as 

^ ^ 1 f,RR ^,1 1 ,oc^ 
i'^ f,R ~ 3 (Af£)2 ' ^ ^ 

where in the last approximate equality we have used the approximate relation ^'^ 



3/, Hi? ■ 

In the regimes C < 1, i-e-, (M^)^ > 1, Eq. ([55]) gives GJ;^™° ~ G/F. In this case m < 1 for sub-horizon 
modes (A: > aH). Thus the deviation from the ACDM model is small, i.e., \F/HF\ < 1 in Eq. ([84]). Hence 
the evolution of matter perturbations is similar to the one in the standard GR case. We recall again that 
this General Relativistic behaviour can be realized even for wbd = because of the presence of a potential 
with a heavy scalar-field mass (M^ ^ k'^ /a?). 

In the regimes ^ > 1, i.e., {Mtf < 1, Eq. ^ gives G^^™^° ~ 4G/3i^. Thus in this case the evolution of 
matter perturbations is different from the one in the GR case because of the appearance of the 4/3 factor. 
If the mass of the Brans-Dicke scalar field is light (Af^ <C fc^/a^), the cosmological effective gravitational 

constant in Brans-Dicke theory is given by G^Jf""° ~ ^ ( 3+2^3° ) 01 ■ Thus in the regime ^ > 1, the f{R) 

theories in the metric formalism behave as the Brans-Dicke theory (with wbd = 0), with a light scalar- field 
mass (M^ < fc^/a^). 



B. Longitudinal gauge (x = 0) 

We shall also derive the approximate equations in the longitudinal gauge {x — 0) foi' sub-horizon modes 
satisfying Eq. ((80|l . We also use the notation a — ^ and tp = — ^I*, which then gives the relation ^> = ^+SF/F 
from Eq. (|75p . In addition to Eq. ((8T|l . we impose the following conditions 

\X\<\HX\, where X = F,F,^,^, (87) 

and 

I'^mXmKm] ■> {H'\B\,H'm,H'\^\} . (88) 

1 LL Cti Lti J 

If the deviation from the ACDM model is not significant, the condition ([87|) is well satisfied. The condition 
([55)1 is also satisfied for sub-horizon modes given in Eq. ([5D|) provided that <I>, ^E* and B are of the same order. 
Under these approximations we obtain, from Eqs. (|67p . ([75|l . (175)1 and ([77)1 . the following relations 

5(„^)+2i7j(?)-f ^$^0, (89) 

'4'^^-^('-^)sp,^, '4^^-^('-±^)sp,^, SF^l^(^)6p^. (90) 

From Eq. ([7^ the term v is of the order of / pm provided that the deviation from the ACDM model is 
not significant. Using Eq. ([90| we find that the ratio 3Hv/ {6p„i/ pm) is of the order of {aH)'^/k'^, which is 
much smaller than unity for sub- horizon modes. This gives (5m ~ Sp„i/pm hi Eq. ()69p . From Eqs. ()89p and 
()90p the matter perturbation in the longitudinal gauge satisfies the following approximate equation 

4^)+2i?<5(?)-g([i|)e)-0. (91) 
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Compared to the comoving gauge the difference appears only in the friction term. Since viable f{R) 
models satisfy the condition \F/HF\ <C 1, Eq. ((84|) reduces to Eq. ([9T|) . We have also checked that in 
uniform density gauge {Spm — 0) the perturbation Sin satisfies the same approximate equation as Eq. (|9ip . 

Before ending this Subsection, we shall introduce a number of parameters which can be useful below. A 
useful parameter is the effective gravitational potential 

$eff = + ^')/2 , (92) 

which characterises the deviation of light rays. This is directly linked with the Integrated Sachs- Wolfe (ISW) 
effect in the CMB [l^, |3l| and weak lensing of distant galaxies [soL [5ll | . From Eq. ([90]) we can approximate 
this parameter by 

'^cff--^^^^^ (93) 



We introduce an anisotropic parameter 



$ - * 2£ 



which behaves as ?7 — > 1 for ^ ^ 1 and 77 — > 2^ for ^ ^ 1. We also define another variable 

I] = q(l+r,/2), (95) 

where q is defined via (fc^/a^)^' = — (1/2) qpmSm'^ . Using the above expressions S can be approximated by 

E-l/i^. (96) 

Note that S is directly linked with $cff- The parameters (S, 77) can be especially important in future survey 
of weak lensing [s^, [5l| . 



C. The appearance of scalarons 

Among the approximations we have used in the previous two subsections, the conditions ([55]) and ([5T|) 
can be violated if an oscillating mode (scalaron) dominates over the matter induced mode discussed above. 
Let us clarify when the oscillating mode becomes important for viable f{R) models satisfying the conditions 
m <C 1 and \F/HF\ ^ 1. For the sub- horizon modes, Eq. (|7^ is approximately given by 

SF + 3HSF +i — + M^\dF~ -Sp^ . (97) 

The solution of this equation is the sum of the matter induced mode SFmd and the oscillatory scalaron mode 
SFosc satisfying 

<5Fosc + SHSFosc + + SFosc = . (98) 

Under the condition {A'P , k'^ / a'^} > H'^ this equation reduces to the form (a^/^JFosc) + uj'^{a^^'^ SFosc) — 0, 
where oj — ^JW^ja^^^rM^ . In the adiabatic regime characterised by ji^/w^l ^ 1 we obtain the following 
WKB solution 

SFosc ^ ca"^/2-^ COS ( / ojdt) , (99) 



where c is a constant. Hence the solution of the perturbation 6R is expressed by 

J,R 1 + 3^ f,RB^V2LJ \J J 
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For viable f{R) models, the scale factor a and the background Ricci scalar evolve as a cx t"^^^ and 
/j(o) ^ ^^^3^2^ during the matter era. Then the amplitude of (5i?osc relative to has a time-dependence 

JliO) °^ (p/a2_,_ ^2)1/4 • ^^"^^ 

Let us consider the models m(r) = C(— r — 1)'' {p > 0) for which the mass M evolves as M cx 
during the matter-dominated epoch. When ^ ^ 1 and ^ » 1 we have |^i?osc|/^^°^ oc t"(3p+i)/2 
|'5^osc|/-R*'°^ oc t~2(p-i-i/3)^ respectively. Hence the amplitude of the oscillating mode decreases faster than 
the background Ricci scalar. Thus if the scalaron is over-produced in the early Universe such that \SR\ > , 
the stability condition (|35p can be violated. This property persists in the radiation-dominated epoch as well 
[2^[29]. Thus in order to ensure the viability of the f{R) theories of gravity in metric formalism, we need to 
esnure that \SR\ is smaller than i?^*'^ at the beginning of the radiation era. This can be achieved by choosing 
the constant c in Eq. ([M|) to be sufficiently small which amounts to a fine tuning for these theories. We note 
that this fine tuning concerns the stability of these theories and is an additional constraint to those usually 
imposed on the parameters of these theories by observations. 

Under the condition that the scalaron mode 5 Rose is negligible relative to the mater- induced mode SR-md, 
one can derive the evolution for the matter perturbation Sm as well as the effective gravitational potential 
$eff- When 5^1 the evolutions of Sm and $off during the matter era are given by 

Sm CX t^/^ , $off = constant . (102) 

Note that the ratio of the matter induced mode relative to the background Ricci scalar evolves as 
I -Rind I /-R'-"-' c>c i^/^ cx Sm- For thc models that satisfy cosmological and local gravity constraints, the 
Universe typically starts from the regime <C 1 and evolved into the regime with ^ ^ 1 during the matter- 
dominated epoch ^28,,^]. When 1, Sm and evolve as 

SmOCt^^-'^', CX ^(^-^)/^ (103) 

For the models m(r) = C(— r — 1)^, we have the time-dependence \5Rind\/R^^^ oc i-2p-i-(\/33-5)/6 ^j^g 
regime ^ ^ 1. This decreases more slowly relative to the ratio \SRosc\/ R^'^'^ oc i~2(p-i-i/3) ^ -j-j^g ggaJaj-on 
mode tends to be unimportant with time. 

In what follows we shall numerically solve the exact perturbation equations in order to check the validity 
of approximations used to reach Eqs. (|84p . (pij) and (|93p. We choose initial conditions such that scalaron 
mode is suppressed relative to the matter induced mode, i.e. \SRlg^\ < \SRl^^^\. We refer the reader to 
Ref. [29] for a comprehensive and detailed study of the scalaron mode. This study also gives the conditions 
under which the scalaron mode dominates over the matter induced mode at the initial stages. 



D. Numerical study of the validity of approximations 



In order to study the dynamics of matter perturbations in the metric formalism we shall introduce the 
following dimensionless variables [ll| 

F f R H ^ 

In terms of these variables, the energy fraction f2„i of the pressureless matter and the effective equation of 
state Wcs are given by 

^m = -^^^^l-Xi~X2~X3, Wcff = -1 - = --^(2x3 - 1) ■ (105) 

The evolution equations for the background dynamics can then be expressed as [ll| 

x[ — — 1 — xs — 3x2 + x^ — xi^a , (106) 
^ ^_^^(^2x3~4~xi), (107) 

TO 

4 = -^-2x3(x3-2), (108) 
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where a prime denotes a derivative with respect to the number of e- folding N = log (a) . For a later use we 
also introduce the variable X4 = aH , which satisfies 



X4 = (X3 - 1)X4: ■ 

The matter epoch corresponds to the critical point 
3™ 1 + 4m 1 + 4m 



Pm ■ {xi,X2,X3) 



1 + m' 2(l + m)2' 2(l + m) 



, WcS 



1 + m 



(109) 



m(7+10m) 
2(l + m)2 ^ 



Since m needs to be much smaller than unity during the matter era, we have Pm ■ {xi,X2,X3) ~ (0, —1/2, 1/2) 
for viable f{R) models. We shall consider a case in which the evolution proceeds from the matter point Pm 
to the de-Sitter point given by: 



Pa - (a;i,a;2,a;3) = (0,-1,2), ^^^ = -1, n,„ = , 
which can easily be shown to lie on the line r = — 2 and is stable for < m < 1. 



(Ill) 



1. Comoving gauge 

In the comoving gauge, the perturbation Eqs. ([78| and ([79|) can be rewritten in terms of the above variables 
thus: 



4^^" + (^3 - Sif - ^(1 - X, -X2- X3)St^ 



— - 3 + 3x1 + 9x2 + 3x3 SF + 3{-2xi + X3 - 1)6F' + 36F" 



6F" + {1-2x1 +X3)SF' + —-X3 
= (1 - XI - X2 - X3)5^Z' - xiSi:^' , 



2x3 



1 — Xl + 3X2 



6F 



(112) 



(113) 



where SF = SF/F. The exact evolution of the matter perturbation can be obtained by solving these equations 
together with the background equations (|106p - (|109p for xi, X2, X3 and X4. Meanwhile, the approximate 
equation ((84|) can be expressed in terms of these variables as 



2^3 



Xi 



2(1 + 30 



Sif - ^(1 - 0:1 - X2 - xa) (^i±|i) 4") , (114) 



where 



(ai/)2 6x3 



(115) 



Let us consider the case in which the condition Af^ ^ k"^ /a^ (i.e., ^ <C 1) is satisfied. Since M needs to be 
large during the matter-dominated epoch to satisfy LGC, this condition holds in viable f{R) models at the 
beginning of the matter era for the modes relevant to large scale structure [H, |2^. Then the term 2x3 /m 
dominates over the term fc^/xl in Eq. (|113p . which gives 5F ~ m5^ under the neglect of scalarons. Hence 
the right hand side of Eq. pi2p can be neglected relative to the left hand side, which means that Eq. I|112p 
reduces to Eq. (|114p . The above argument shows that, in the regime C ^ Ij Eq. (|114p can be valid even for 
super-Hubble modes as long as the contribution of the scalaron is unimportant. In this regime the matter 



perturbations evolve as in the case of standard GR, i.e. 5\ 



[v) 
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FIG. 1: The evolution of the matter perturbation 5™' in the comoving gauge for the model m(r) — (—r — l)"* with the 
mode k/aoHo — 10. Initial conditions were chosen to be xi — 0, X2 — —0.5000, xs — 0.5001, (5™' = 10~^, (5m' = 10"'^, 
5F = 8.0 X 10"^^ 5F' = and k/aiHi = 4.1 at the redshift z = 28.9. The solid curve is obtained by solving the 
exact equations l|112p and pi3[) numerically, whereas the dotted one is obtained by solving the approximate equation 

(Uni). 



The perturbations can enter the regime <C fc^/a^ il^-^ C ^ 1) before reaching the present epoch, 
depending on the mode k and on the evolution of M (28l l29j . In models m(r) = (— r — 1)^ this occurs for 
the modes k/aoHo > 3.5, where the subscript represent present values. In the case k/aoHo — 300, the 
redshift at k/a = M corresponds to = 4.83. Since is always larger than in the past because of the 
requirement m ^ 1, the modes are inside the Hubble radius [k^ /o?H^ > 1) after the perturbations enter 
the regime JVP < k'^ /a? . Hence the approximation we used to reach Eq. (|114[) is valid in this regime. In the 

regime < fc^/a^ the term (fc^/a;|)(5F in Eq. (|113p . balances the term {\ — xi — X2 — x^)5m\ which gives 
rise to an additional contribution on the right hand side of Eq. (jll2p . This then leads to the approximate 
equation (|114p with ^ ^ 1, which has a growing-mode solution (5^ oc i(v^-i)/6^ 

In Fig. [1] we plot the evolution of 5m^ for the model m(r) = (— r — 1)'^ with the mode k/aoHo = 10. 
Initial conditions are chosen so that the scalaron mode does not dominate over the matter-induced mode. 
In this case the transition from the regime Af^ > k'^/a^ to the region < k'^/a^ occurs at the redshift 
Zk = 1.62. We find that the approximate equation pi4p shows an excellent agreement with the results 
obtained by numerically solving the exact equations (|112p and (|113p . The argument also holds for modes 
that are initially outside the Hubble radius. Thus the approximate equation ()84p is reliable to estimate the 
growth of matter perturbations and the resulting matter power spectrum, provided that the scalaron does 
not dominate in the early Universe. 

2. Longitudinal gauge 

In the longitudinal gauge the combination of Eqs. ((73)) - (|77p leads to the following perturbation equations 

+ (^2 - ^xi + X3^ $' + (3x2 + 3X3)$ = lx2SF - Qx3 + 1^ SF' - ^SF" , (116) 

SF" + (X3 + 2)5F' + fl^ + 3X2 + —) SF = fexa + 2x3 - $ - (3xi + 2)$', (117) 

V 3 X4 m J \ 3 X4 / 



(2 + 3x2 -X3 + 2x3/m + P/xl)6F + (2 + xi + xs)SF' + SF" + (xi - 6x2 - 2x3)$ + (4xi + 2)$' 

1 - Xi - X2 - X3 
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(118) 

^ - 2$' + (2 - + 5F' + (1 + xi)6F , (119) 
where we have used 5* = $ + dF. The effective potential defined in Eq. ([92|l is given by 

^eS = ^ + ^5F. (120) 

In order to understand the evolution of perturbations at the initial stages of the matter era, let us consider 
the regime ^ ^ 1 without assuming the sub- horizon condition k/{aH) ^ 1. We have in mind viable f{R) 
models with vanishingly small values of m deep inside the matter epoch. Equation (|117|) then becomes 



5F ~ -2m 



$-2m$'. (121) 



Note that under sub-horizon approximation we have 5F ~ —2mk^^/3{aH)^, which agrees with Eq. (|90|) . 
Using Eq. p2ip we find that the right hand side of Eq. (|116D can be neglected relative to the left hand 
side, thus giving the solution $ = constant (together with a decaying mode proportional to t^^/"^). From 
Eqs. (IMl) and (fTT^ we obtain $cff ^ * and 

2^.2 
3[aH )^i2m 

4^)' ^ 4^) . (123) 

Equation p22)) agrees with the expression (j93[) obtained under the sub-horizon approximation (fc/a 3> i?). 

Since $cff is a constant, the matter perturbation can be seen from Eq. (|122[) to evolve as Sm^ oc a. This is 
consistent with the approximate equation (j9ip . i.e., 

-f xs^)' - ^(1 -x,-x,- X3)ii|^4^ ^ , (124) 

which has the growing mode solution Sm^ — Sm^ cx a in the regime ^ ^ 1. 

We may ask why the above method reproduces the result derived under the sub-horizon approximation, 
without employing the approximation k/a 3> H. In the regime ^ 1 the perturbation SF is suppressed 
relative to $ as given in Eq. (|12ip . This allows us to neglect the right hand side of Eq. pi6p . giving a 
constant $. This mimics the situation in General Relativity where 6F = and $ = constant together with 
Eq. (|122p . resulting in Sm^ on a. Moreover, from Eq. (|119[) . the quantity B — Hv -I- \E' is well approximated 
by B ~ 5^/3 = constant. Hence the right hand side of Eq. ([57)1 can be neglected even without assuming 
the sub-horizon approximation. Thus using the relation p22p we can obtain Eq. (j9ip in the regime ^ -C 1 
without assuming k/a 3> H. The above approximation corresponds to the limit of large M (M^ ^ k^/a^), 
which gives rise to the evolution of perturbations close to the case of General Relativity. In General Relativity 
{6F = and F — 0), one has the exact equation (|122p from Eqs. ([75]) and ([7^ . Thus the perturbations in 
the large M case (^ ^ 1) mimic those in General Relativity, apart from the fact that the scalaron is present 
in the former but not in the latter. 

When ^ ^ 1 one has k^ ja} ^ ^ which means that the sub-horizon type approximation we used 
in the subsection B holds well in this regime. This situation is similar to the case of the comoving gauge. 
For the modes that start from the regime Af ^ ^ fc^/a^ and enter the regime Af^ ^ k^ jd} before the end of 
the matter era, the evolution of perturbations changes from the standard general relativistic form (jl02p to 
the non-standard form (|103p . 

In Fig. [2] we plot the evolution of b'^m and <I>cff in the model m(r) = {—r— 1)"^ for the mode k = aoHo that 
lies outside the Hubble radius at the start of integration (z — 28.9). Together with numerically integrating 
Eqs. (|116p - (|118p . we also solve the approximate equation (|124p with <i>off derived by ([55)1 . From Fig. [5] we 
find that the approximate equations agree well with the exact numerical results, even if the mode is initially 
slightly outside the Hubble radius. We caution, however, that for large-scale modes far outside the Hubble 
radius the scalaron can be important. In fact we have numerically checked that the oscillating mode appears 
for such super-Hubble modes unless the coefficient of the scalaron in Eq. (|100p is fine-tuned to be small. In 
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Fig. [2] the growth of the gravitational potential is not seen in the region z < Zk, since the transition redshift 
is small (zk ~ 0.36). It can, however, be observed if we consider modes on smaller scales. 
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FIG. 2: The evolution of 5™' and $cft in the longitudinal gauge for the model m(r) — (— r — 1)'' with a mode 
k = aoHo- We numerically integrated Eqs. (|116|) . (|117|l and (|118|l with initial conditions <I>cff — 10"^, O^ff = 
and k/aiHi — 0.36 and with 5^' and J^' satisfying Eqs. (|122p and p23p . Initial conditions for the background 
quantities were chosen to be the same as in Fig. [T] We also plot 5m ' and <^a« obtained by solving the approximate 
equations (|124|l and (|93p . The approximation is valid even when the mode is initially outside the Hubble radius. 

In summary, for viable f{R) models that satisfy the cosmological and local gravity constraints, the ap- 
proximate Eqs. (|124p and are trustable even for the modes outside the Hubble radius, as long as the 
scalaron is suppressed relative to the matter-induced mode. 

E. Constraints on the model m(r) = C{—r — 1)^ 

Finally we consider the current and future constraints on models m(r) — C{—r — 1)^ with < C < 1. At 
the background level, compatibility with the SNIa observations could result in the divergence of the equation 
of state of dark energy [12, |23| ■ Interestingly the redshift at which such a divergence may occur could be of 
order unity. However the current SNIa observations are not yet sufficiently accurate to rule out such cases. 
Some constraints on the model parameters can be obtained from the present equation of state of dark energy, 
but even models with p — 1.5 and C — 1 are allowed [29]. Thus the background does not provide strong 
constraints on the model parameters. However this situation can change in the future when higher-redshift 
data will become available from the observations of SNIa and Gamma Ray bursts. 

There are a number of additional observational constraints on the growth rate of matter perturbations. At 
the redshift z ~ 3, McDonald et al. [s^] obtained the constraint (5^/5™ = 1.46 ±0.49 from the measurement 
of the matter power spectrum from the Lyman-a forests. Taking into account the more recent data reported 
by Viel and Haehnelt [s^ in the redshift range 2 < z < 4, the maximum value of the growth rate allowed 
by the current observations is given by [53 | 

s = 5'J5„,<l.b. (125) 



The current data still have large error bars and some data even allow the parameter range — 1 < s < 2 [53| . 
However it is expected that in future observations the growth rate will be constrained more severely. From 
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Eqs. HnH) and (fTM)) we have s = 1 for > fc^/a^ and s = (\/33 - l)/4 = 1.186 for < P/a"^. In 
Fig. [3]we plot the evolution of the growth rate for models m(r) = (— r — 1)"^ for a number of different values 
of k. The increase of s from unity implies that the perturbations enter the regime ^ k"^ /a^ . For smaller 
scale modes this transition occurs earlier, which leads to the larger maximum value of s. The growth rate 
begins to decrease once the Universe enters the late-time accelerated epoch. As estimated analytically, the 
growth rate is bounded by s < 1.186. Hence the current observational constraint p25p is still too weak to 
place constraints on m(r) — C{—r — lY models. 
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FIG. 3: The evolution of the growth rate s = S'^/Sm with respect to the redshift z in the model m{r) — (— r — 1)'^ 
with four different values of k. Initial conditions were chosen as in Fig. [S] The transition redshift zt defined as the 
redshift where k/a — M, becomes larger for smaller scales. After the matter perturbation enters the region z < Zk 
the growth rate begins to increase toward the value s = 1.186, but it starts to decrease once the Universe enters the 
stage of accelerated expansion. 



However, these models exhibit peculiar features in the matter power spectrum. This is a consequence 
of the fact that there is a transition redshift Zk at which the growth rate begins to change from s = 1 to 
s = 1.186. For the modes relevant to galaxy clusters {k/aoHo = 0(10^)), this transition typically occurs 
during the matter-dominated epoch (see Fig. [3]). Since the time at z = depends upon the modes k 
(tk oc fc~3/(3p+i)-j^ ^]-^jg leads to the change in the slope of the matter power spectrum. The difference between 
the slopes of the matter power spectrum determined from galaxy surveys and the CMB spectrum, on the 
scales k/aoHo = 0(10^), is given by 



V 33 — 5 , , 

An.^-^. (126) 

This analytic result agrees well with numerical results except for models with p 3> 1 [29j . Observationally no 
significant differences have so far been found between the two power spectra. If we take the bound An < 0.05 
we obtain the constraint p > 5. To place further constraints on models, a likelihood analysis is required 
which employs the data from both the galaxy power spectrum and the CMB. Such observational constraints, 
including the data from the Lyman-a forest and Gamma Ray bursts, are under consideration 55]. 

Numerically, we find that models m(r) = (— r — 1)^ have constraints on the present value of the deviation 
parameter given by m(z = 0) ^ 10^^. Thus even though m is constrained to be very small during the matter 
era, a notable deviation from the ACDM model can occur around the present epoch. 

Finally the ISW effect in the CMB power spectrum is important on large scale modes with fc/oo-ffo — C(l)- 
As can be seen from Fig. [21 even models with p = 3 and C = 1 do not give rise to a significant amplification 
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of the gravitational potential. The models with p > 2 are consistent with the low multipoles in the CMB 
data [23 • Thus this effect does not generally provide stronger additional constraints. 



V. DENSITY PERTURBATIONS IN THE PALATINI FORMALISM 



In this section we discuss the evolution of density perturbations in the Palatini formalism and the resulting 
observational consequences. The background equations are given by 

FR-2f = -p™ , 
where the Ricci scalar R satisfies the following relations 

3 



(127) 
(128) 

(129) 
(130) 

The matter perturbations satisfy Eqs. (|63p and (|64p as in the case of the metric formalism. Other pertur- 
bation equations are [47 1 



R = 6{2H^ + 
3Hp, 



H) + — [ F + 3HF - — 
' F \ 2F 



R=-- 



RFi 



fc2 
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' 2F 



Ha — if 



_1_ 3F' 
2F \ 2F 



3HF a 



2F 



SF- (n+^UF-Fa + p^v 



X + Hx-a~<f^j{SF-Fx) 



(131) 

(132) 
(133) 
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F 



3F2 

'f^ 



2F F'^ 



3F . 

a 

2F 



- R + — \ SF + isH \ SF + 3SF 



(134) 
(135) 



Note that we corrected several typos found in Ref. [47]. Given the non-dynamical nature of Eq. ()135|1 . it 
is clear that the scalaron mode does not exist in the Palatini case. This is associated with the fact that 
the Palatini formalism corresponds to generalised Brans-Dicke theory (Eq. PH)) ) with cjbd ~ —3/2. The 
perturbation SF is directly determined by the matter perturbation Spm, as 

F^R Sprn 



SF - 

F 1-m 

where m is defined in Eq. (|14p . 

Below we shall consider these perturbations in various gauges. 



(136) 



A. Comoving gauge 

In the comoving gauge (v = 0) one has a = and k = Sm^ . Then from Eq. (|134p we find 
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1 

2F 




3iJ - — Uf + 3(5F + 5p„ 



(137) 



Combining this with Eq. (|136p we obtain an exact second-order equation for 5m\ which can be numericaUy 
solved (See Appendix A for details). On the other hand, since we are mostly interested in the evolution of 
modes on sub-horizon scales, it makes sense to consider the approximate equations similar to those considered 
in the metric case. 

Using a sub-horizon type approximation, such that only the terms containing k"^ / a? and 5pm are considered 
on the right hand side of Eq. (|137p . together with Eq. (|136p . we obtain the following approximate perturbation 
equation 

'^'Z' + f + ^] C'^ - (l + Sl::^ - , (138) 



where ^ is defined in Eq. 

Alternatively we may study the case in which the deviation from the ACDM model is small, i.e., 

|m|<l, (139) 

as required from LGC The derivative terms such as |to'| and \m"\ are also assumed to be much smaller 
than unity. Using the fact that from Eqs. (|128p and p36p the perturbation SF in this case is of the order of 

mSm\ Eq. ()138p can be obtained under condition (I139P without using the sub-horizon approximation. Thus, 
if the deviation from the ACDM model is small, the approximate equation (|138p is valid even for the modes 
outside the Hubble radius. This situation is similar to the case in the metric formalism. In fact we have 
confirmed this property by numerically solving the exact equation p37p and comparing it with the solutions 
of the approximate equation (|138p . 

One can estimate the order of the term R6F on the r.h.s. of Eq. (I134p by using Eq. (|136p . i.e., RSF = 

mSpm/i^ ~ '^)- This gives rise to the contribution of the order of {p„i/2F)m6rn' in the third term of 
Eq. (|138p . which is negligible under the condition (|139p . As long as we neglect this contribution, we should 
approximate ^/(l — m) ~ ^ in the third term of Eq. (|138p . In the following, we implicitly assume this when 
we write the term (1 — m) in the denominator. 

In the limit ^ = -^rj^rn ^ 1, Eq. (|138p agrees with Eq. (|84p in the metric formalism. However a significant 
difference appears in the regime ^ S> 1. In that case there is a strong amplification of the matter perturbation 
in the Palatini case because of the growth of the term ^ in Eq. I|138p . We shall estimate this growth rate for 
a number of concrete models in Subsection C below. 



B. Longitudinal gauge 



We next consider the Longitudinal gauge (x = 0), and as in the metric case we use the notation a — ^ 
and Lp = Under the sub-horizon type approximation used in the comoving case above, the evolution 

equation reduces to Eq. (|89p obtained in the metric case. Using Eqs. (|13ip and (|133p together with Eq. (|136p . 
these approximations also give 

-$^-—1 + -^ Up™, c^-—[l--^\5pm. (140) 

2t \ 1 — my zi* \ 1 — m/ 

Hence the matter perturbations satisfy the following approximate equation 

^L^) + 2i/<j(^) - (l + 4^^ - . (141) 

The effective gravitational potential ^eff defined in Eq. ((92)l satisfies 



- 2fc2 F ™ ' 



(142) 
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which is the same as in the metric case. Similarly the parameters rj and S defined in Eqs. ([M)) and ([M]) in 
this case become 

V^ , S^l. (143) 

1 — m — ^ 

We note that while the expression for rj is diff'erent from that in the metric case, S remains the same. 

The above approximate equations (|14ip and p42p are valid under the conditions (|139p and C 1 even 
without the sub-horizon approximation. In fact the argument is similar to the metric case in which Eqs. (j9ip 
and ([Mjl reduce to the corresponding GR equations for ^ <C 1. 

We also note that in regimes ^ ^ 1 the perturbation modes are inside the Hubble radius, which shows 
that the sub-horizon approximation above is still valid. Thus, as long as the condition (|139p is satisfied, 
we can safely use Eqs. (|141l) and (|142p even for super-Hubble modes. Also since in the Palatini case the 
perturbation 5R is sourced by the matter induced mode only, we do not need to worry about the dominance 
of the scalaron oscillations for super-Hubble modes, unlike in the metric case. 

In Ref. [4?! the equation for matter perturbations was derived in the uniform density gauge {5pm = 0). 
This is an exact equation as in the comoving gauge (see also Ref. t56i|)- We shall show in Appendix B that 
similar to other gauges Eq. (|139p is sufficient to reduce the exact equation to the approximate one. We shall 
also show that the approximate equation is the same as the one in the longitudinal gauge (|14ip . 

Since the evolution of matter perturbations hardly depends on the gauge chosen, we shall in what follows 
denote the matter perturbations simply by 5^- 



C. Analytic estimate for the growth of perturbations 



As we mentioned above the evolution of perturbations in the regime ^ ^ 1 is similar to the standard GR 
case where 5m oc i^^'^, s = 5'^/ 6m = and <I> cx constant. In this subsection we shall estimate the growth 
rate of perturbations after the system enters the regime C > 1- We shall consider models with |m| <C 1, for 
consistency with LGC P5|) . 

During the matter era in which the Ricci scalar evolves as i? oc the parameter ^ is given by ^ = 
±ma/mkak, where the subscript "fc" denotes the values when the system crosses ^ = 1. Here we note that 
the plus sign corresponds to a positive m, whereas the negative sign to a negative m. As we showed in 
Sec. II, the latter case can be allowed unlike the metric case. Under the condition |m| <C 1, the matter 
perturbation (|14ip satisfies the following equation 



ISL - - f 1 ± —e^-^" 



0. 



•^m ' 2 ' 2 

Lets us consider the case in which the evolution of the parameter m is given by 

where p is a constant. The values of p in several different /(i?) models are given by 
• (i) f{R) = aR^+"' - A (constant m): p^O, 

p^l- P for R:$> Rc, 

WR^Y" . p=2n+l for R->Rc, 



(144) 



(145) 



• (ii) fiR)^R^XRc 
. (iii) /(i?) = i?-Ai?c(^^75^ 



(iv) f{R) =R~XR, 



1 



1 + ^ 



p = 2n+l for i? > i?c 



With the above choice of m, Eq. (|144p reduces to 



0^ 



2 2 



lie 



{ap+l){N-N^) 



0. 



(146) 



For the positive sign in Eq. (|146p . i.e., for m > 0, the solution of Eq. (I146P can be written in terms of a 
linear combination of Bessel functions and Y^, : 



(147) 
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where ci , C2 are constants and 

ygg(3p+l)(Ar^7V,)/2 5 

X = , V ^ . 148 

For the negative sign in Eq. (|146p . i.e., for to < 0, the solution of Eq. (|146|) is given by 

5^ = e-(^-^'=)/4 [ci J, (x) + C2Y, {x)] , (149) 

where x and v are as given in (jl48p . 

In the foUowing we shall discuss the positive and negative to cases in turn. 



1. m > 



As an example, we consider the constant to models {p — 0). In this case, the Bessel function J^{ix) has 
a growing mode solution J^/2{ix) oc I^/2{x) cc e^/^/x for a; 3> 1 (here l5/2{x) is a modified Bessel function 
with X = VEe'-^-^"'^/^). Then, in the regimes ^ S> 1, the evolution of the matter perturbations and its 
growth rate are given by 

5™ oc exp(V6e(^-^'-)/2) , ^ = = ^e(^-^'=)/^ (150) 

where we have used -\/6e(^"^'=)/^ ^ {N — Nk)/2. Thus the growth rate s of the matter perturbations 
increases very rapidly. Also from Eq. (|142p . in the regimes ^ ^ 1, the effective gravitational potential grows 
(double) exponentially as 

$efr oc exp(e^(^-^^)/2) , (151) 

which leads to a strong and observable ISW effect. 

Similarly, in models with p ^ 0, one can estimate the evolution of perturbations in the regime ^ ^ 1: 

<5„, oc c&eff OC exp ( 3^^^ j , . = V6^(3,+i)(^„^,)/2 _ ^^^^^ 

This shows that for models with p > the growth rate increases faster than in the constant to models. 
When p < —1/3 the above instability can be avoided, but in that case to increases towards the past. Thus 
unless the present value of m is negligibly small, the condition |to| <C 1 required for LGC can be violated 
during the matter era. Thus these models are hardly distinguishable from the ACDM model in the present 
Universe. In what follows we shall concentrate on the positive p case. 

We note that for these positive to models the violent growth of Sm and <I>cff comes to an end after the 
Universe enters the accelerated stage, since ^ and flm begin to decrease. 



2. m < 

When m is negative, the Bessel functions in the regime |^| 3> 1 behave as Jiy{x) ^ ^Ij-Kx cos[x — {2v + 
l)7r/4] and Y^(x) ^ ^^jux sin[a;— (2(^+1 )7r/4] respectively. Thus the solution of Eq. (|149|) in this asymptotic 
region becomes 

<5™~Ce-(3p+2)(A^-^'=)/4cos(a; + 0), (|e| » 1) , (153) 
where C and 9 are constants. Using this solution, we obtain 

^ -i(3p + 2)5,„-^Ce3p(^-^'')/4sin(x + 0), (154) 
s ~ -h3p + 2)- ^^^^ xta.n{x + e), (155) 
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which are also vaHd in the regimes |^| ^ 1. When p > 0, Sm exhibits damped oscillations whereas 
increases in time with the oscillations. The averaged value of the growth rate s is given by s = — (3p + 2)/4, 
but it shows a divergence every time x changes by tt. 

If the Universe crosses the critical point |<^| = 1 around the end of the matter era, it does not necessarily 
reach the regime |f | ^ 1. In such cases one can not fully use the above approximate solutions. We shall 
show later that, in some cases, the Universe can enter the accelerated stage without oscillations in S„i up to 
the present epoch. The oscillations in S^, can be seen as we choose larger values of |rn| and k. The frequency 
of oscillations tends to grow for larger values of p. The models that enter the regimes |^| 3> 1 generally have 
difficulty in being consistent with observations, since they lead to largely negative values of s as given by 
Eq. USSl). 



3. Constraints on \m\ from the requirement \^\ < 1 

The f{R) models can be consistent with observations if the Universe does not enter the regime |^| > 1 
until the end of the matter-dominated epoch. One can estimate the ratio of the comoving Hubble radius 
(aH)^^ during the matter era to its present value thus: 



aoHo ( a 



1/2 



where c = 1 in the absence of the dark energy dominated epoch. The presence of a dark energy era leads to 
a change in the value of c. Numerically this factor is around c— 1.7-1.9. Using the relation R ~ 2>H^ that 
holds during the matter era for |m| <C 1, we find that |^| crosses 1 at a critical redshift 

If Zc is smaller than the order of unity, the Universe does not enter the regime |^| > 1 during the matter 
dominated epoch. This gives the following constraint: 

|m(z)| < (^^) , for z>0{l). (158) 

The matter power spectrum, in the linear regime, has been observed for the scales 0.02/i Mpc^^ ^ ^ 
0.2/iMpc""'^. Non-linear effects can be important for smaller scales with k > 0.2ft, Mpc~^. Taking the value 
k ~ 0.2ft Mpc~^ ~ GOOoo-ffoj below which linear perturbation theory is valid, we obtain the constraint 
|m(z)| ^ 3 X 10^^ during the matter era. 

Of course this is a rough estimate and the actual constraints on m{z) depend upon the particular models 
considered. For example, even if |^| crosses 1 during the matter era, the models can be consistent with 
observations provided that |^| does not grow rapidly after the crossing. Whether or not ^ reaches the regime 
1^1 3> 1 depends on the particular models chosen. Hence to place constraints on m, we need a detailed 
analysis for each particular model. In the next subsection we shall provide a numerical investigation of 
a number of f{R) models presented above and place constraints on present values of m as well as model 
parameters. 



D. Constraints on model parameters 



In this subsection we shall employ the information provided by the growth of the matter density pertur- 
bations to place constraints on the parameters of the f{R) models presented in subsection C above. We do 
this by numerically solving the exact perturbation (|165p together with the background equations (|127p and 
(|128p . See Appendix C for equations written in a form convenient for numerical integration. 



1. Constant m models: f{R) = 0/?^+™ - A 

Compared to other models considered here, the growth of |^| is rather mild in the constant m models, 
being of the form |^| oc a = . Thus, in order for |^| to grow from 0.1 to 10, one would require an increase 
in the number of e- foldings by 4.6. 
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We shall first consider the positive m case. In the left panel of Fig. |3]we plot the evolution of growth rate 
s = S[^/5m for the mode k = 600aoi?o for several values of m. For m = 3 x 10~^ we numerically obtain 
Zc ~ 11, denoted by a black dot in Fig. [H This almost agrees with the analytical estimate (|157p which gives 
Zc ~ 10. In the regimes ^ 1 the evolution of matter perturbations is given by 5™ = 5^ oc a, which results 
in s ~ 1. The growth rate s begins to move away from unity as f becomes of order 0.1, and then continues 



to grow before the Universe enters the stage of accelerated expansion. For this model we find s„ 



2.06 



and ^max ~ 3.13, which shows that the model does not enter the regime ^ ^ 1 where the evolution of 
perturbations is described by Eqs. pSOp and (|15ip . 

For a model with m = 1.5 x 10~^, the critical redshift occurs at around Zc ~ 5 with s ~ 1.4. The maximum 
value of the growth rate is Smax ~ 1.57, which corresponds to the marginal case satisfying the criterion (|125p . 
For a model with m = 2.0 x 10~^, the evolution of perturbations is not much different from the general 
relativistic case. 
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FIG. 4: The evolution of perturbations for the model: f{R) ~ aR^'^"^ — A with positive values of m. In the left 
panel we show the growth rate s = S'^/Sm versus the redshift z for the mode k/aoHo = 600 with three different 
values of m. The black dots represent the points at which ^ crosses 1. The right panel depicts the evolution of s for 
m — 1.5 X 10~^ with three different values of k. 



To show the variation of the growth rate as a function of scale, we depict in the right panel of Fig. [3] the 
evolution of s for the model m — 1.5 x 10^^ for three different values of k. As can be seen, the maximum 
value of the growth rate s decreases as k is decreased (i.e. the scales become larger). In particular, for 
the mode k = lOOaoi/o (corresponding to fc = 0.33ft. Mpc~^), the evolution of perturbations exhibits no 
difference compared to the corresponding evolution in the general relativistic case. Hence the matter power 
spectrum is enhanced on small scales (fc — 0.1h-0.2h Mpc~^), while the spectrum remains similar to the 
standard general relativistic case on larger scales (fc = 0.02/i-0.04ft Mpc~^). This results in different spectral 
indices on different scales. Placing more precise constraints on m, would require performing a likelihood 
analysis using the data from the matter power spectrum. However, in order to obtain an order of magnitude 
for the maximum value of m, it is sufficient to use the criterion (|125p for the mode fc = 600ao_ffo. For the 
constant m models we find the constraint to be m ^ 10~^. 

When m is negative, the growth rate s decreases unlike the positive m case. In the left panel of Fig. [S] we 
plot the evolution of s for three different negative values of m for the mode fc/oo-ffo = 600. As can be seen s 
tends to decrease more rapidly with increasing \m\. If m = —2.0 x 10~^ the present value of s becomes very 
small (s < —1). As we see in the right panel of Fig. [5l when m — —2.0 x 10"*^, there is a significant fall in 
the values of s between k/aoHa = 300 and k/aoHo — 600. This can lead to large differences in the spectral 
indices of the matter power spectrum for small and large scale modes. From the above argument |m| should 
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be smaller than the order of 10~^, which has an upper bound similar to the positive m case. 

When TO = —2.0 x 10~^ the Universe crosses the point |^| = 1 at the redshift Zc ^ 7.4, but the increase of 
1^1 for z < Zc is mild. Moreover the quantity |^| begins to decrease after the Universe enters the accelerated 
stage. Numerically we obtain the value ^ ~ —0.77 at present (z — 0). Thus the system does not reach the 
regime |^| 3> 1, and hence not a single period of oscillation occurs by the present epoch. However, for larger 
values of |m|, we have numerically checked that the oscillations of b^a indeed occur. 
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FIG. 5: The evolution of perturbations for the model: /(i?) = Q.R — A for negative m. In the left panel we show 
the growth rate s = 5^/5™ versus the redshift z for the mode k/aoHo — 600 with three different values of m. The 
black dots represent the points at which the quantity |^| crosses 1. The right panel depicts the evolution of s for 
m = —2.0 X 10~^ with three different values of k. 



We also recall that the growth of effective gravitational potential $off leads to an ISW effect in the CMB 
spectrum on large scales {k/aoHo r^a, few). However, when \m\ ~ 10~^, $off does not grow for these modes. 
As a result the ISW effect does not provide stronger constraints to than those provided by the matter power 
spectrum. 



2. fiR) =R- \Rc{R/Rcf 

The observational constraints on the parameters of this model were studied in Ref. [i^ . [Note that Rc is 
not very much different from the present value of the cosmological Ricci scalar Rq.] Here, we shall obtain 
constraints on the parameter m which for this model is given by 

A/3(l-/3)(i?/i?c)^-i 

and make a comparison between our results. The late-time de-Sitter point {R ~ Ri) is obtained from the 
constraint equation FR — 2/ = 0, to give [Ri/ R^Y^^^ = A(2 — 0). Thus at this de-Sitter point the variable 
TO satisfies 

m(i?i)=/3/2. (160) 
For /3 < 1, the parameter m in the regime i? S> i?c is given by 

TO ~ A/3(l - I3){R/R,f'^ oc , (161) 
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which decreases towards the past. 

If /3 (< 1) is of the order of unity, the quantity m is too large to satisfy the requirement (|158p for the mode 
k = 600ao-ffo during the matter era (recall that from Eq. (|160p the present value of m is of the order of 
/3). This is basically associated with the fact that, in the regimes R ^ Rc, the model gives a linear relation 
between m and r [m — C{—r — 1)]. Thus we need the condition \f3\ ^ 1 in order to be compatible with the 
criterion (|158p . 

To determine the changes in the behaviour of this model as a function of /3, we considered three distinct 
values of /3 and calculated the corresponding growth rate s and the parameter m in each case. Our results 
are summarised in Fig. [6l The left hand panel shows the evolution of the growth rate s for A = 1 and 
k = 600ao-ffo with the three different values of (3. For /3 = 1.5 x 10~^ the present value of the parameter m 
is around m(zo) ~ 6.7 x 10"^, which is close to the value of m at the de-Sitter point (m(i?i) = 7.5 x 10~^). 
We also find that the parameter f crosses 1 at a redshift Zc ~ 3 with m{zc) ^ 1.2 x 10^*^. 




FIG. 6: The evolution of perturbations for the model }{R) = R — XRc{R/Rc)^ with positive /3 and A = 1. The left 
hand panel depicts s — S'^/Sm versus the redshift z for the mode k/aoHo = 600 with three different values of p. 
The right hand panel shows the evolution of m with respect to z for k/aoHo = 600. From the requirement (|125|l we 
obtain the constraint /3 < 8.2 x 10""". 

Furthermore, we find that the growth rate s is larger for these models than in the case of constant m 
models. This is due to the fact that ^ in this case evolves faster, as ^ cx ^2(4/3-/3) r^^ie maximum growth 
rate reached for (3 = 1.5 x 10^** corresponds to Smax ~ 1-88 with ^ ~ 4. As expected, models with smaller 
values of /3 possess growth rates which are more compatible with observational constraints. Employing the 
criterion (|125p for the mode k — GOOgoHq, we find the constraint f3 < 8.2 x 10~^. This is slightly larger than 
the constraint /3 < 3.0 x 10~^ obtained in [3^ from the likelihood analysis of the SDSS data [53|- In the left 
panel of Fig. [6] we also consider this case in order to find the corresponding evolution of s. The maximum 
value of the growth rate in this case is found to be Smax ^ 1.095, which indicates that the constraint (|125p is 
rather weak. Nevertheless, the criterion ()125p is certainly sufficient in order to extract the order of magnitude 
of the bound on /3. 

In the right panel of Fig. [6] we plot the evolution of the parameter m for the case with A = 1 and 
k = 600aoi?o, for three values of (3. As can be seen m increases from the past to the present. Using the 
criterion (|125p we obtain the bound m{z = 0) < 3.5 x 10~^. If we adopt the severer criterion s < 1.1, the 
constraint becomes rn(z = 0) < 1.3 x lO^'^. Thus the deviation from the ACDM model is constrained to be 
small (m(z = 0) < IQ-^). 

We also studied the effects of changing the parameter A in the action on the bounds on (3. We considered 
two cases with A = 10 and A = 100. We found that these changes in A have negligible effects on the 
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constraints imposed on (3 and m{z = 0), compared to that obtained from the case with A = 1. The reason 
for this lack of sensitivity is that a change in parameter A is compensated by corresponding changes to the 
values of Ri , oq and Hq . 

When /3 < the parameter m is negative from Eq. (|16ip . In the left panel of Fig. [7] we plot the evolution 
of s for three different values of f3 with k/aoHo — 600. We find that the present values of s become smaller 
than —1 for \(3\ > 1.2 x 10^*, in which case \m{z = 0)| is smaller than the order of 5.3 x 10~^ (see the right 
panel of Fig.[71). Thus if we use the criterion s{z = 0) ^ —1 for the validity of the models, the upper bounds 
of |/3| and \m{z = 0)| are similar to those in the positive /? case. 
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FIG. 7: The evolution of perturbations for the model f{R) = R ~ XRc{R/ Rc)^ with negative /3 and A = 1. The left 
hand panel depicts s = S'^/Sm versus the redshift z for the mode k/aoHo — 600 with three different values of The 
right hand panel shows the evolution of the quantity m with respect to z for k/aoHo = 600. If we use the criterion 
s{z = 0) > —1, we obtain the constraint /3 > —1.2 x 10"*. 



3. f{R) = R- \Rc\l- {1 + R" /Rl)-"\ 

Finally we consider the above model (where n > 0) recently discussed by Starobinksy [28j . The parameter 
m for this model is given by 



2nAa;(l + x2)-"-2[(2n + l)a;2 - 1] 
1 - 2nAa;(I + 



^^r^i;.. ^^X": ^ -^ere .^R/R.. (162) 



and the de-Sitter point at i? = i?i corresponds to 

A = — iEi(l + a^i) — where xi = Ri/Rc. (163) 

2[(1 + j;2)"+i - 1 - (n+ 1)j:2 ' i i/ c k j 

Once we fix the value of A, x\ is known accordingly. In the regime R^ Rc the parameter m behaves as 

m ~ 2n(2n + l)A(i?c/i?)^"+^ oc t2(2n+i) _ (164) 

Because of the presence of a larger power of {Rc/R) in the expression for m in this case, m decreases more 
rapidly towards the past compared to the model f{R) ~ R — \Rc{R/ RcY^ discussed above. For the mode 
k — 600ao-ffoj the bound (|158p implies that m has to be smaller than the order of 10^^-10"^ by the end of 
the matter-dominated epoch, in order to ensure that the model does not enter the regime with ^ > 1. 
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In Fig. [5] we plot, for mode k = 600aoi?0: the evolutfon of s and m for A = 2.5 with three different values 
of n. When n = 3.07 the critical redshift is given by Zc ^ 1.05 with m ^ 1 x 10~^. The rapid increase of 
s occurs in the regime ^ > 1, after which the growth rate reaches a maximum value Smax ~ 2. The present 
value of m is found to be m = 4.5 x 10~^, which is an order of magnitude larger than its corresponding 
value at ^ = 1. Using the criterion (|125p . we obtain the constraints n > 3.23 and m{z = 0) < 2.9 x 10^"* 
for A = 2.5. The present value of m in this model is one order of magnitude larger than the corresponding 
values in the constant m models as well as f{R) = i? — \Rc{R/ Rc)^ model. 




FIG. 8: The evolution of perturbations for the model f{R) ^ R- \Rc[l - (l + R^/Rc) "] with A = 2.5. The left 
hand panel depicts s = 5'm/Sm versus the redshift z for the mode k/aoHo = 600 with three different values of n. The 
right hand panel shows the evolution of the quantity m with respect to z for k/aoHo ~ 600. 

We also find that in contrast to the model f{R) ~ R—XRc{R/Rc)^ the constraints on n for the Starobinsky 
model are sensitive to the values of the parameter A. For larger values of A the constraints on n is weaker. For 
example, for A = 10 and A = 50 we find the corresponding constraints on n imposed by (|125p to be n > 1.74 
{m(z = 0) ~ 1.4 X 10^*) and n > 1.09 (m(z = 0) 1.1 x 10^*) respectively. This can be understood in the 
following way. When A is increased, we obtain a larger ratio Ri/Rc from Eq. (|163p . which also leads to a 
larger ratio R/Rc in the past. Then from Eq. (|164p a smaller value of n is sufficient to realize the condition 
|m| ^ 1. It can also be seen from the form of the action that the values of Rc can also affect the constraints 
on n. We find that for small A values, Rc has a small effect on the constraint, whereas for large values of A 
the affect of changing Rc is negligible. 

From Eq. (I164|) we find that m can be negative for — 1/2 < n < (and A > 0) in the regime R 3> Rc- 
When n is close to 0, the models are close to the model f{R) — R ~ XRc{R/ Rc)^ discussed above. We 
find that s{z ~ 0) is larger than —1 for \n\ < 9.3 x 10^^, in which case we have \m{z = 0)| < 4.5 x 10^^. 
When n is close to —1/2, Eq. (|164p seems to suggest that the models should be close to the constant m 
models. However, care needs to be taken in this case since m changes sign from negative to positive at 
(R/Rc)'^ ~ l/{2n + 1) in the deep matter dominated epoch. As a result, for n close to —1/2, we numerically 
find that the growth rate s shows a rapid growth for (R/Rc)^ < l/(2n + 1). Thus, in the limit n — > —1/2, 
the models do not behave as constant m models and they are excluded observationally. 

We have also analysed the model f{R) = R — XRc (jiJji^yL_^_i {n > 0) of Hu & Sawicki [27| and have found 
the constraints on the parameters n and m{z = 0) to be n > 3.33 and m(z = 0) < 2.15 x 10~^ respectively 
for A = 2.5 and k/aoHo = 600, which are similar to the constraints derived above. 

In summary, the present values of m are constrained to be m{z = 0) < 10^"* from the bound p25p in both 
Starobinsky and Hu & Sawicki models. 
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VI. CONCLUSIONS 

We have made a detailed study of the evolution of density perturbations in f{R) gravity theories in both 
metric and Palatini formalisms and employed them to study the viability of models in each case. To study the 
viability of concrete models we considered three sets of constraints, provided by the background cosmological 
evolution, local gravity experiments and the evolution of matter density perturbations respectively. 

We began by considering the cosmological and local gravity constraints. For models satisfying these con- 
straints, we proceeded to study the additional constraints provided by the evolution of density perturbations 
to further constrain the model parameters as well as their deviation from the ACDM model. 

The f{R) theories in the metric formalism are equivalent to generalised Brans-Dicke theories with a 
scalar-field potential V{(j)) and Brans-Dicke parameter cjbd = 0. The presence of the field potential, allows 
the construction of f{R) models that satisfy the local gravity constraints under the use of a chameleon 
mechanism. We find that for typical models of the forms ([?7)l and to satisfy the cosmological and local 
gravity constraints, the parameter m is required to be much smaller than unity during the radiation and 
matter eras but can grow to values of order of 0.1 in the accelerated epoch. Models in the metric formalism 
also suffer from an additional fine tuning due to the presence of scalaron oscillating modes (which is absent 
in the Palatini case). Finally, to be stable these theories require J^br to be positive. 

On the other hand, the f{R) theories in the Palatini formalism correspond to generalised Brans-Dicke 
theories with a scalar-field potential V{(f>) and Brans-Dicke parameter lobb = —3/2. This makes these 
theories special in the sense that the oscillating scalar degree of freedom (scalaron) is absent in these theories 
and therefore the corresponding fine tuning to the metric case does not exist. Also unlike the case of the 
metric formalism, there is in this case no notion of field mass M that determines an interaction length 
mediated by a fifth force. Thus the LGC for these theories need to be analyzed separately in contrast to 
theories with lubd ^ —3/2. The main condition required in this case in order to satisfy the LGC is that |r7i| 
is smaller than the order of unity. Moreover the requirement for the cosmological viability in the Palatini 
formalism is not severe compared to the metric case. Thus in contrast to the metric case, to satisfy the 
cosmological and local gravity constraints, we do not require vanishingly small values of m during radiation 
and matter dominated epochs and furthermore j ^rr can be negative in this case. As a result, even models 
of the type /(i?) = R- ^2("+i)/i?" with n > can be allowed at the background level unlike the metric 
case. 

We then studied the constraints provided by the evolution of density perturbations in each case. In the case 
of the metric formalism we derived the equations for matter perturbations under sub-horizon approximations 
in several different gauges. In regimes ^ fc^/a^ (i-S-, C 1)j found the approximate perturbation 
equations to be valid even without using sub- horizon approximations, provided that the scalaron mode is not 
dominant relative to the matter-induced mode. This is a consequence of the fact that, when ^ /c^/a^, 
the evolution of perturbations mimics that in General Relativity. After the Universe enters the regime 
<C fc^/a^, the modes are inside the Hubble radius due to the fact that the condition M"^ ^ Ris required 
for the compatibility with LGC. Thus, for the models that satisfy LGC, as long as the scalarons do not 
dominate over the matter-induced mode, approximate perturbation equations are valid even for the modes 
that initially lie outside the Hubble radius. In the Palatini case the approximate equations are even more 
reliable because of the absence of scalarons. 

In the metric formalism, most viable f{R) models take the form m(r) — C{~r — 1)^" {p > 1) in the 
regimes where the Ricci scalar is larger than the order of the present cosmological value. In these models, 
the modes relevant to the observed matter power spectrum correspond to the regimes ^ fc^/a^ with 
the growth rate s — S'„^/6m = 1 at the beginning of the matter era. These models typically enter the regime 
<gC k"^ /a? during the matter era in which the growth rate of matter perturbations is given by s = 1.186. 
If we use the present observational bound s < 1.5, we do not obtain strong constraints on these models. 
However, since the transition time at k/ a = M depends upon the mode fc, there is a difference in the spectral 
indices between the matter power spectrum and the CMB spectrum [see Eq. (|126l) ]. If we take the bound 
An < 0.05, the models with p > 5 are allowed. The present value of the parameter m is constrained to be 
m(z = 0) ^ 10~^. Thus, while m needs to be negligibly small during the radiation and matter eras, one can 
have appreciable deviation from the ACDM model around the present epoch. 

In the Palatini formalism the approximate matter perturbation equations are valid even for super-Hubble 
modes, for models satisfying LGC (|m| <C 1). If m is positive, there is a strong amplification of 5^ in the 
regime ^ ^ 1, whereas if m < the matter perturbation exhibits a damped oscillation for |^| ^ 1. When 
the quantity m evolves as m oc during the matter era, we have analytically estimated the growth rate 
s in both positive and negative m cases [see Eqs. (jl52p and (jl55l) . respectively]. From the requirement 
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that the Universe does not enter the regimes |^| > 1 during the matter era, we obtain the constraint 
|m(z)| < (aoHo/k)^ for z > 0(1). While this is a good criterion to avoid non-standard evolution of matter 
perturbations, one needs to carry out a more detailed analysis to place constraints on the quantity m for 
each f{R) model. When m is positive, we have obtained the constraint m ^ 10~^ by considering the modes 
k relevant to the matter power spectrum. We also studied the evolution of perturbations for the models 
f{R) = R-XR,{R/Rcf and f{R) = i?- Ai?c[l - (1 + i?V^c)""]- For these models we found the constraints 
m{z = 0) ^ 10~^ and m{z = 0) ^ 10~*, respectively, from the requirement s ^ 1.5. Thus, unlike the metric 
case, the deviation from the ACDM model at the present epoch is small even when m grows from the matter 
era to the accelerated epoch. This situation does not change for negative values of m. 

In summary, for viable models in the metric formalism, the quantity m is constrained to be very much 
smaller than the order of unity during the matter era from LGC, but it can grow to the order of 0.1 around the 
present epoch. In the Palatini formalism, LGC and background cosmological constraints do not place strong 
bounds on m (only requiring |m| ^ 10~^), but the density perturbations can provide stringent constraints: 
|m| < 10^^-10^^. Thus in the Palatini case the f{R) theories are hardly distinguishable from the ACDM 
model even at the present epoch. This follows from a peculiar evolution of the matter perturbations in the 
Palatini case, in the regime |^| > 1, that exhibits rapid growth (when m > 0) or damped oscillations (when 
TO < 0). 

While the constraints obtained here are sufficient to give the orders of magnitude of the allowed model pa- 
rameters, it will be of interest to obtain more precise constraints by using recent and upcoming observational 
data including large scale structure, CMB, Supernova la, gamma ray bursts and weak lensing. 
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Appendix A: The equation for matter perturbations in the comoving gauge in the Palatini formalism 



In this appendix we present the exact matter perturbation equation in the comoving gauge in the Palatini 
formalism. As in the metric case this equation needs to be solved simultaneously with the background 
equations (|127p - (|128[) . Unlike the metric case, however, it is not easy in this case to find dimensionless 
variables in terms of which both sets of equations close. As a result we proceed to integrate the equations 
directly. Using the background equations and ignoring the radiation, the perturbation equation (|137p can 
be written as 



Pl<5L^)"+P2<5L^)'+P35L"^=0 

where the coefficients Pi, P2, P3 are given by 
3J 



(165) 



Pi = 
P2 = 
Pi = 
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Wn''' F 




(168) 
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where J is defined as 



All the terms in the coefficients Pi , P2 , ^3 can be expressed in terms of the scale factor a (or equivalently 
N), which thus allows Eq. (|165p to close and be readily integrated numerically. 



Appendix B: The equation for matter perturbations in the uniform density gauge in the Palatini 

formalism 



In the uniform density gauge {5pm = 0) we have 6F = = 5R from Eq. (|136p and v^^^ — a, k = 
SHv^^l + Kv^^'^ from Eqs. ^ and dH]). Substituting these relations for Eq. pM)) we obtain 



3 H 



F_ 
2F 



3F 3HF 3F2 
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F2 



where we have used the relation 



R = 6(2ij2 + h) + -{f + SHF - 



F 



El 

2F 



2F a2 



=0, 



Then the matter perturbation, Sm^ — '3Hv''-^\ satisfies the following equation of motion 



where 



Cl 



C2 = 



2H 



1 + F/2HF 
1 + F/2HF 
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m / 2HF 2iJ2F2 2H^F 
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(170) 



(171) 



(172) 



(173) 



H 2H H I F 



F H H 



\ H F 

llF \ ~ 2IP ^ 6a2iJ2 I " Jpjpp 



(174) 



This agrees with the result derived in Refs. [43, |5y| . 

Let us obtain the approximate equation for matter perturbations under the approximation (|139p . Taking 
note that the term \F /HF\ is of the order of |m|, the coefficients ci and C2 are given by 



ci^2H, C2 = 



H 2H F 



i?3 i?2 QHF {aHY 



(175) 



When we estimate the first two terms in the square bracket of C2, we use the following approximate 
relations 



2FiJ ~ -p™ , 2FH ~ 3Hp,n , 
which come from Eqs. ([T?f)) and From Eq. wc have 



F ■ 



F -RF 



,R 



(176) 



(177) 



Using these relations, wc find that the matter perturbation satisfies the following approximate equation of 
motion 



This is the same equation as the one in the longitudinal gauge (|14ip . 



(178) 
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Appendix C: Equations convenient for numerical simulations in The Palatini formalism 



In this appendix wc shall present the equations convenient for numerical simulations. From Eqs. (|127p . 
(fT28l) and (fnoj) we obtain 



= „ — , where C 



2p„, +FR-f 



6FC 

Introducing a dimensionless quantity 



3 FR{FR-2f) 
2 F{F,rR - F) 



2 



(179) 



FR~ f , , 



we get the differential equation for y [23]: 

y' ^y(l-y)[3 + C{H)] , (181) 



where C{R) is defined in Eq. 
The following relations also hold 



FR-f 2y 



FR-2} l~y 



(182) 



^ = 1 - y ■ (183) 

Specifying the value of y, the initial Ricci scalar R is determined by Eq. (|182p . Solving Eq. (|18ip . we obtain 
?/, i?, Fl and from Eqs. (|182p . (|180p and (|183p . The effective equation of state of dark energy is given by 

F C FR / ^ 

^,, = -y+__ + _^__^. (184) 

As long as the deviation from the ACDM model is small (|m| ^ 1), we have WeS — —yi- 
The perturbation equations (|14ip and (|142p are given by 

C + ^(1 - 3w.fi)SL - ^C(i ~y){i + Srn^O, (185) 

3 faH^ ^ 



'i'cfr ^ "2 j C(l-y)'5m. (186) 

Although we solve exact perturbation equations, the above perturbation equations are found to be very 
accurate. 
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